Dataset2f_100_100sp

Same as Dataset 2f but sample size of 100 and number of taxa to 100. All different interaction matrices

1. Load Library

---

library(tidyverse)

── Attaching core tidyverse packages ────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     ── Conflicts ──────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(igraph)

Attaching package: ‘igraph’

The following objects are masked from ‘package:lubridate’:

    %--%, union

The following objects are masked from ‘package:dplyr’:

    as_data_frame, groups, union

The following objects are masked from ‘package:purrr’:

    compose, simplify

The following object is masked from ‘package:tidyr’:

    crossing

The following object is masked from ‘package:tibble’:

    as_data_frame

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

library(NBZIMM)

Attaching package: ‘NBZIMM’

The following object is masked from ‘package:stringr’:

    fixed

library(SpiecEasi)

Attaching package: ‘SpiecEasi’

The following object is masked from ‘package:igraph’:

    make_graph

library(LIMON)
library(here)

here() starts at /Users/bealab/Desktop/Long_networks

library(lme4)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:SpiecEasi’:

    tril, triu

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(Matrix)
library(tscount)
library(patchwork)
library(MASS)

Attaching package: ‘MASS’

The following object is masked from ‘package:patchwork’:

    area

The following object is masked from ‘package:SpiecEasi’:

    fitdistr

The following object is masked from ‘package:dplyr’:

    select

library(matrixcalc)

Attaching package: ‘matrixcalc’

The following object is masked from ‘package:igraph’:

    %s%

library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

library(devtools)

Loading required package: usethis
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(miaSim)

Loading required package: TreeSummarizedExperiment
Loading required package: SingleCellExperiment
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count


Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads,
    colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds,
    colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2,
    colTabulates, colVarDiffs, colVars, colWeightedMads,
    colWeightedMeans, colWeightedMedians, colWeightedSds,
    colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    rowWeightedSds, rowWeightedVars

Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following object is masked from ‘package:gridExtra’:

    combine

The following objects are masked from ‘package:igraph’:

    normalize, path, union

The following objects are masked from ‘package:lubridate’:

    intersect, setdiff, union

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position,
    rank, rbind, Reduce, rownames, sapply, setdiff, table, tapply,
    union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:Matrix’:

    expand, unname

The following objects are masked from ‘package:lubridate’:

    second, second<-

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following object is masked from ‘package:utils’:

    findMatches

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following object is masked from ‘package:lubridate’:

    %within%

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians

Loading required package: Biostrings
Loading required package: XVector

Attaching package: ‘XVector’

The following object is masked from ‘package:purrr’:

    compact


Attaching package: ‘Biostrings’

The following object is masked from ‘package:base’:

    strsplit

library(reshape2)

Attaching package: ‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

library(corpcor)

Attaching package: ‘corpcor’

The following object is masked from ‘package:matrixcalc’:

    is.positive.definite

2. Simulate Counts

Simulate a different starting covariance matrix and counts per subject. ie each one has their own seed

# Initalize an empty list
####################################################
count_tables1 <- list()

# The Loop
####################################################

# Set the indices
for (i in 1:100) {
  
  # Set the seed for each subject
  set.seed(12345 + i)  
  
  # Create Starting covariance matrix
  n <- 100^2  
  lower_values <- runif(n / 2, min = -0.7, max = -0.2)  
  upper_values <- runif(n / 2, min = 0.2, max = 0.7)
  user_interactions <- c(lower_values, upper_values)
  user_interactions <- sample(user_interactions)
  user_A  <- randomA(
    n_species = 100,
    interactions = user_interactions,
    diagonal = -1.0,
    connectance = 0.5,
    scale_off_diagonal = 0.2,
    symmetric = TRUE)
  
  # Heatmap of the covariance matrix
  color_palette <- colorRampPalette(c("red", "white", "blue"))(20)
  heatmap(
    user_A,
    Colv = NA,
    Rowv = NA,
    col = color_palette,
    zlim = c(-1, 1), 
    main = paste("Subject", i)
  )
  
  # Run the simulation
  tse_glv <- simulateGLV(n_species = 100,
                         A = user_A,
                         t_start = 0, 
                         t_store = 4,
                         stochastic = FALSE,
                         norm = FALSE,
                         error_variance = 0.1)

  # Get the count table
  sim_data <- tse_glv@assays@data@listData[["counts"]]
  
  # Store the count table in the list
  count_tables1[[i]] <- t(sim_data)
}

# Merge together
####################################################

# Combine the count tables
combined_count_table <- do.call(rbind, count_tables1)

# Add Rownames
rownames(combined_count_table) <- paste0("Sbj", rep(1:100, 
                            each = nrow(count_tables1[[1]])), "_Time", 1:4)

3. Create fake metadata

Make Metadata and merge with the count data

# Set seed
set.seed(12345)
# Df 1 is Metadata
########################################################
meta_data <-  expand.grid(Time = 1:4,ID = 1:100)
rownames(meta_data) <- rownames(combined_count_table)



# Df 2 is Metadata merged with Counts
########################################################
#Round off and increase
combined_count_table <- as.data.frame(combined_count_table + abs(min(combined_count_table)))
combined_count_table <- (combined_count_table)*10
meta_counts <- base::merge(meta_data, combined_count_table, by ="row.names", 
                           all = TRUE)
meta_counts <- column_to_rownames(meta_counts, "Row.names")

Plot the data to look at the differences among subjects

Graphs to Check

# Individual Species Plots
########################################################
# Pivot to long data
count_long <- tidyr::pivot_longer(meta_counts, cols = starts_with("sp"), names_to = "Species")

# Plot the data
count_long %>%
  ggplot(aes(x = Time, y = value, colour = as.factor(ID),
             group = as.factor(ID), linetype = as.factor(ID))) +
  geom_line() + 
  geom_point() +
  geom_jitter() +
  ylab("Count") +
  labs(linetype = "ID", color = "ID") +
  facet_wrap(~ Species) +  # Create a panel for each species
  theme(legend.position = "none") +
  ggtitle("Time Series Data")

# Plot only two
########################################################
# Pivot to long data
metacounts_filt <- meta_counts[,c(1:6)]
count_long <- tidyr::pivot_longer(metacounts_filt, cols = starts_with("sp"), names_to = "Species")

# Plot the data
count_long %>%
  ggplot(aes(x = Time, y = value, colour = as.factor(ID),
             group = as.factor(ID), linetype = as.factor(ID))) +
  geom_line() + 
  geom_point() +
  geom_jitter() +
  ylab("Count") +
  labs(linetype = "ID", color = "ID") +
  facet_wrap(~ Species) +  # Create a panel for each species
  theme_linedraw() +
  theme(axis.text.x = element_text(family = "arial",color = "black", size = 14),
        axis.text.y = element_text(family = "arial",color = "black", size = 14),
        axis.title.x = element_text(family = "arial",color = "black", size = 14),
        axis.title.y = element_text(family = "arial",color = "black", size = 14),
        legend.position = "none",
        strip.text = element_text(face = "bold", family = "arial", color = "white", size = 15))

# Distribution of counts
########################################################
hist(as.matrix(combined_count_table), 
     breaks = 100, main = "Distribution of GLV Data", xlab = "Counts")

4. Save the data

Save the Counts

# Make more like counts
meta_counts[,3:102] <- round(meta_counts[,3:102])
write.csv(meta_counts, here("Data","GLV_SimData", "Dataset_2f", "GLV_100_100sp.csv"))

# Save just timepoint 2
meta_counts_filt <- meta_counts %>% filter(Time == 2)
write.csv(meta_counts_filt, here("Data","GLV_SimData", "Dataset_2f", "GLV_100_100sp_T2.csv"))